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Abstract 

We propose a method for multi-scale hybrid simulations of molecular dynamics (MD) and com- 
putational fluid dynamics (CFD). In the method, usual lattice-mesh based simulations are applied 
for CFD level, but each lattice is associated with a small MD cell which generates a "local stress" 
according to a "local flow field" given from CFD instead of using any constitutive functions at 
CFD level. We carried out the hybrid simulations for some elemental flow problems of simple 
Lennard-Jones liquids and compared the results with those obtained by usual CFDs with a New- 
tonian constitutive relation in order to examine the validity of our hybrid simulation method. It 
is demonstrated that our hybrid simulations successfully reproduced the correct flow behavior ob- 
tained from usual CFDs as far as the mesh size Ax and the time-step At of CFD are not too large 
comparing to the system size ^md and the sampling duration £md of MD simulations performed 
at each time step of CFDs. Otherwise, simulations are affected by large fluctuations due to poor 
statistical averages taken in the MD part. Properties of the fluctuations are analyzed in detail. 
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I. INTRODUCTION 



Hydrodynamics of complex fluids are of particular importance in various science and en- 
gineering fields, such as fluid mechanics, soft matter science, mechanical engineering, chem- 
ical engineering, and so on. Because of the complicated couplings between internal degree 
of freedoms of complex fluids and their flow behavior, conventional treatments, based on 
usual assumptions such as non-slip boundary conditions and linear Newtonian constitutive 
relations, are often invalid. Striking examples can be seen in systems such as colloidal disper- 
sions, polymeric liquids, granular matters, and liquid crystals. Those systems are known to 
exhibit peculiar flow behaviors, e.g., shear thinning or thickening, viscoelasticity, jamming, 
flow induced phase transition, etc. 

Although there exits huge accumulation of experimental and theoretical studies on the 
rheology of complex fluids, performing computational fluid dynamics (CFD) simulations are 
not yet common for complex fluids since reliable constitutive equations are often unknown 
for those systems. On the other hand, there exists a different problem also for microscopic 
approaches such as molecular dynamics (MD) simulations, while constitutive equations are 
no more necessary in this case. The characteristic time and length scales of complex fluids 
easily become several orders larger than those of microscopic scales. Therefore, most hydro- 
dynamic problems of complex fluids are yet out of reach of microscopic MD simulations. To 
overcome those serious limitations mentioned above, we aim to develop a new multi-scale 
method which is for performing hybrid simulations of MD and CFD valid for complex fluids 
without any constitutive equations. 

Various methods for hybrid simulations of MD and CFD has already been proposed by 
several researchers. Most of those methods are based on "domain decomposition" for which 
MD simulations are applied only around the points of interest, i.e., in the vicinity of defects, 
boundaries, interfaces, where details of molecular motions are important, while the remaining 
regions are treated only by CFD.^jMiMiLSiSJii Exchange of information between MD and 
CFD is performed in a coupling regions where each system is subjected to some constrains 
to take the consistency of the two systems. This kind of hybrid method is expected to be 
useful specially for problems including interfaces, such as adhesion, friction, anchoring of 
crystal liquids, stick-slip motions, etc. 

In order to apply hybrid methods of MD and CFD to hydrodynamics of complex fluids, 
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a different type of approach is probably needed. Our strategy for this is straightforward. 
We try to develop a multi-scale hybrid method based on the local equilibrium assumption. 
Here CFD is used as a fluid solver, while MD simulations are used only to generate local 
properties, such as constitutive relations of the fluid under consideration, by performing local 
statistical sampling in a consistent matter. The numerical algorithm is rather simple. We 
perform usual lattice-mesh based CFD simulations at an upper level, but each mesh-node is 
associated with a small lower level MD cell which passes a "local stress" to CFD according 
to a "local flow field" given from CFD to MD instead of using any constitutive functions at 
CFD level. MD simulations thus have to be performed at all node points and at every time 
steps of CFD. 

One might think that the simulations would be be much faster if we construct tabular 
database of the constitutive relations by performing MD simulations in advance under many 
different simulation parameters and refer the table from CFD. The "tabular approach" works 
much effective for simple fluids for which the constitutive relations depend only on a few 
parameters, such as the density, temperature, and shear rate. In the case of complex fluids, 
however, the number of parameters to be considered can be huge depending on the local 
quantities to be considered. In the case of charged systems for example, the local stress 
depends also on the local compositions and chemical potentials of ions and the local electric 
field, etc. Although we used only simple Lennard- Jones liquid in the present study, we adopt 
the local sampling strategy rather than the tabular contraction strategy to be more general. 
The main purpose of the present study is to examine the validity of our multi-scale hybrid 
model by performing some simple demonstrations of the method. Efficiency and drawback 
of both strategies will be considered in future for more specific problems. An idea similar in 
spirit to the present method was also put forward earlier by W. Ren and W. E.— 

The hybrid simulation method is described in Sec. II, and some demonstrative results 
for one- and two-dimensional flows of simple Lennard- Jones liquids are shown in Sec. III. 
A special attention is put on the efficiency and the reliability of our hybrid method there. 
The simulation results obtained by our multi-scale hybrid method are compared with those 
of normal CFDs with a Newtonian constitutive relation. The validity of our method is 
discussed in Sec. IV, and a summary is given in Sec. V. 
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II. HYBRID MODEL 



Incompressible flows for isotropic materials are described by the following equations, 

dv " 0, (1) 



dv a dv a ldP a p 

where x a is the Cartesian coordinate system, t the time, v a the velocity, p the density, P a p 
the stress tensor, and g a the external force per unit mass. Here and after the subscripts a, 
(3, and 7 represent the index in Cartesian coordinates, i.e. {a, /3, 7} = {x, y, z}, and the 
summation convention is used. The stress tensor P a p is written in the form, 



a/3 



-pS a/3 + T a/3 , (3) 



where p is the pressure and 8 a p is the Kronecker delta. Here we assumed that the diagonal 
component of the stress tensor is isotropic. The off-diagonal stress tensor is symmetric T a p = 
Tp a and traceless T aa =0.— In order to solve the above equations, one needs a constitutive 
relation for the stress tensor T a p. In our hybrid method, instead of using any explicit formulas 
such as the Newtonian constitutive relation, T a p is computed directly by MD simulations. 



A. CFD Scheme 



We use a lattice-mesh based finite volume method with a staggered arrangement for 
vector and scalar quantity— See Fig. [TJ The control volume for a vector quantity is a unit 
square surrounded by dashed lines and that for a scalar quantity is a unit square surrounded 
by solid lines. Eqs (pQ) and are discretized by integrating the quantities on each control 
volume. As for numerical time integrations, we use the fourth order Runge-Kutta method, 
where a single physical time step At is divided into four sub-steps. More concretely, the 
time evolution of a quantity 0, which is to be determined by the equation d(f>/dt=f(t, 0), is 
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FIG. 1: Staggered arrangement of vector quantity, the velocity v, and scalar quantity, the pressure 
p and density p, on a lattice-mesh grid. 



written as 



(,** 



+ ^/(*n,<A (4a) 



At/fo^,^), ( 4c ) 

/(t„,0™) + 2/(t„ +l ,0; +i )+ 

» 2 

2/(^+1/2, 0n+l/2) + /(Wl) 0n+l)] ■ 

(4d) 



Time evolution of the fluid velocity i? is computed by the above set of equations. On the 
other hand, the pressure p is determined so that the fluid velocity satisfies the incompressible 
condition ([1]) at each sub-step. The procedure at each sub-step is written as 

p = p + ijj, (5a) 
v = v — r'Vtp, (5b) 

= -V«, (5c) 

T 

where p is the pressure obtained at the previous sub-step, v is the velocity obtained by 
solving equation (jlj) at the present sub-step, and r is the time increment of the sub-step. 
The remaining three components of the tensor T a p are to be computed directly by MD 
simulations. The detail of the method is described in the next subsection. Note that the 
calculations of T a p is carried out at each sub-step of equation (jlj). 
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(a) Space dDcompasLt-ion 



(b) Time evolution 



FIG. 2: Schematic diagram for the hybrid scheme, (a) CFD simulations are performed in a reference 
coordinate (x,y,z), while MD simulations are performed in a rotated coordinate (x',y',z') so that 
the diagonal components of E' a p become all zero with the procedure described in Sec II B. The CFD 
system is discretized into cubic subsystems whose side length is Ax. Each subsystem is associate 
with a MD cell, whose side length is Imd, with Lees-Edward periodic boundary condition under 
shear deformation, (b) A schematic time evolution of our multi-scale method. CFD simulation 
proceeds with a time step of At, MD simulation is carried out for a lapse of time £md only to 
sample local stress T' a a at each node point and time step of CFD. 

B. Computation of Local Stress by MD 

We compute the local stresses by MD simulations according to the local strain rates, 
rather than the local flow velocities themselves, computed at the CFD level. A schematic 
diagram of the method is depicted in Fig. [2j At the CFD level, the local strain rate tensor 
E a p is denned as 



2 \dx/3 dx ay 

where the incompressible condition, E aa =0, is to be satisfied. We can now define a rotation 
matrix with which the strain rate tensor E a p is transformed to 

/ 

E' = QEQ T = 



pi 



E 'xy 


E 'xz 


\ 











E 'yz 








E 'zy 





/ 



(7) 



where the diagonal components all vanish. This transformation makes performing MD sim- 
ulations much easier with the usual Lees-Edwards periodic boundary condition for simple 
shear flows under the assumption that each off-diagonal the component of the local stress 
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tensors depends only on the corresponding component of the local strain rate tensors, re- 
spectively. The off-diagonal stress tensor T' a g is computed according to E' a/3 and then pass to 
CFD after transforming back to the original coordinates, T a p. For one- or two-dimensional 
flows [d/dz=0 and f^=0], and E' are expressed as 

(cos 8 sin 9 \ 
(8) 
— sin 9 cos 9 I 

E' xy = E' yx = —E xx sin 29 + E xy cos 20, (9) 

where 

e =\^ {-ty < io) 

Non-equilibrium MD simulations for simple shear flows in the rotated Cartesian coor- 
dinates are performed in many MD cells according to the local strain rate E n s defined at 
each lattice node of the CFD. The number of particles in each MD cell is 256 if not men- 
tioned. Once a local stress tensor P' a a is obtained at the MD level, the local stress at each 
lattice node P a p in the original coordinate system is obtained by combining the pressure p 
obtained a priori by CFD and a tensor T' a a obtained by subtracting the isotropic normal 
stress components from P' a p as 

p = e T [- p i + T'}& = -pi + e T T'e, (11) 

where I is the unit tensor. For one- or two dimensional flows, we can use T' =T' =0 and 

xx y y 

rpl _ r pl _ TDl 

xy yx xy 

In the non-equilibrium MD simulations, we use the Lees-Edwards sheared periodic 
boundary condition to a cubic MD. The temperature is kept at a constant by using a 
thermostat.— 1 ^ The stress P' af5 is averaged in steady states after transient behavior van- 
ished. 



III. NUMERICAL COMPUTATION 



We have carried out the hybrid simulations for one- and two-dimensional flows of a simple 
liquid composed of the Lennard- Jones (LJ) particles interacting via the potential 



v LJ (r) 



Ae 



r . 



r 



(12) 
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In the present simulations, the potential is truncated at r=r c and sifted to zero at the dis- 
tance for computational efficiency. We considered only the cases where the temperature 
T and the fluid density p are uniform and constant over the CFD systems and the ex- 
ternal force is neglected, g a =0. The reduced temperature T*=Tk/e and reduced density 
p*=pa 3 /m, where k is the Boltzmann constant and m is the mass of a single LJ particle, are 
fixed at T*=1.0 and p*=0.8 in the simulations. Here and after, non-dimensional quantities 
normalized by the energy and length parameters of the Lennard- Jones potential, e and a, 
are denoted by the superscript "*" . 

In the following, At and Ax represent the time-step and the mesh size of CFD calcu- 
lations, and £md and /md represent the sampling time and the side length of a MD cell, 
respectively. The two parameters At/t MD and Ax/Imd represent the efficiency of our hybrid 
simulations. We have carried out hybrid simulations with several different values of the 
parameters and compared the results with those obtained by usual CFDs. In the present 
simulations we fixed t MD =0.005 and / MD =6.84, while At and Ax are changed as listed in 
Tabled 

A. Pressure-driven channel flows 

The Lennard- Jones liquid with r*=2.5 is contained in channel composed of two parallel 
plates located at x\=±L/2 and subjected to a pressure gradient in y-direction. We performed 
one- and two-dimensional simulations for this pressure- driven channel flows. The pressure 
gradient is set as Ap/(pU 2 / X)=1.25, where Ap is the pressure difference over a distance L, 
and U is a characteristic flow velocity. Non-slip boundary condition is applied on the two 
plates. 

The results of one-dimensional simulations are shown in Figs. |3] and HI A symmetric 
condition is used at x=0, and the computational domain [— L/2,0] is divided into eight slits. 
Parameters used in the simulations are listed as SP I — III in Table [B In the corresponding 
CFD simulation, viscosity is set as i]*=2.0, which is for the LJ liquid with r* = 2.5 at T*=1.0 
and p*=0.8. The Reynold number defined as pUL/rj is fixed at 40. It is clearly shown in 
Fig. [3] that the results obtained by the present hybrid simulations well agree with those of 
usual CFDs for the case At/tMD < 2 and Ax//md < 2. When the parameters become large, 
instantaneous velocity profiles tend to fluctuate as seen in Fig. HI It should be noted that 
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FIG. 3: The velocity profiles obtained by one dimensional computations for the pressure-driven 
channel flow. Simulation Parameters are summarized as SP I in Table [I] for (a) and SP II in Table 
UJfor (b). The solid lines show the results of usual CFD simulation, and dotted lines and square 
symbols show the results of the present hybrid simulation. 

the fluctuations can be removed almost perfectly by taking time averages. This means that 
the mean values of the fluctuation is almost zero, i.e., the fluctuation might be removed also 
by applying some filtering, etc. We will discuss on this point later. 

The results of two-dimensional simulations are shown in Fig. [5j The computational do- 
main is now [-L/2,L/2] x [0,L/2] and divided into 16x8 uniform lattices. Non-slip boundary 
condition is used at x=±L/2. At y=0 and L/2, periodic condition is used for the velocity, 
and the pressure is set as p(x, L/2)=p(x, 0) — 0.5Ap. Here we assumed Stokes flow, i.e., the 
second term of the left-hand-side of Equation <^ is dropped in the computation. It is seen 
that, at each time step, the velocity fluctuations are much smaller than the pressure fluc- 
tuations. This is clearly due to the incompressible condition to be imposed to the velocity. 
The velocity also fluctuates immediately after solving Eq. (j2j), however the incompressible 
condition Eq. (TTJ) tends to adjust it. The pressure fluctuations can be removed also by 
taking time averages. 

B. Two-dimensional cavity flows 

Lennard- Jones liquid with r* = 2 1 / 6 is contained in a square box whose side length is 
L. At t—0, the upper wall starts to move from left to right at a velocity v w =U. Non-slip 



9 




FIG. 4: The velocity profiles obtained by one dimensional computations for the pressure-driven 
channel flow. Simulation parameters are summarized as SP III in Table HI Time evolutions of 
the velocity profile after an application of pressure gradient in y-direction at i=0. Instantaneous 
profiles at t/At = 200, 500, and 1500 are shown in (a), while the velocity profiles are time averaged 
over t/At = [0,300], [300,600], and [1200,1500] in (b). The squares show results of the hybrid 
simulations, and the dotted lines present the corresponding CFD results for comparison. 




(a) instanUaTLBOLia (b) i irnoiivcrragecL 



FIG. 5: The steady-state flow profiles of the pressure-driven channel flow obtained by a two- 
dimensional calculation. Simulation parameters are summarized as SP IV in Table HI 

boundary condition is applied at each wall; v x =U and v y =0 at y=L, and v x =v y =0 at other 
walls. At left- and right- upper corners, v x =U and v y =0 is applied. The results of the 
hybrid simulations are shown in Figs. M and [7J The computational domain is divided into 
32x32 uniform lattices. Values of parameters used in the present simulations are listed as 
SP V-VII in Table HI The Reynolds number is defined as pUL/r], and the viscosity of the 
corresponding LJ fluid is rf—1.1. 

Fig. [6] shows the steady-state velocity profiles time-averaged over t/At— [950,1000]. It is 
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(a) R==5g (b) H==235 (□] Rc=941 



FIG. 6: The steady-state velocity profile for the cavity flow. Simulation parameters are summarized 
as SP V in Table UJ for (a), SP VI for (b), and SP VII for (c). The velocity profiles are time averaged 
over t/At= [950, 1000] 

clear that our hybrid method can successfully reproduce the characteristic flow properties 
of cavity flows with different Reynolds numbers. Fig. [7] shows time evolutions of velocity 
profiles for the case of Re=980 after a sudden application of upper-wall sliding at t — 0. 
Here, the results obtained by hybrid simulations are compared with those of usual CFD 
simulations. It is seen commonly that a small vortex first appears at the upper-right corner 
is moving gradually toward the center of the box with increasing the size of the vortex as time 
passes. The agreements between hybrid simulations and CFDs are very well. Our hybrid 
method is confirmed to reproduce successfully the time-evolution while large fluctuations 
are seen in the instantaneous velocity profiles. 

IV. DISCUSSION 

As mentioned above, the ratios At/t MD and Ax/Imd measure the efficiency of our hybrid 
simulations. Larger the ratios, simulations are more efficient, however, the statistical fluc- 
tuations also become large. For example, in a case of At/tuD=Ax/luD=^, computational 
efficiency is, roughly speaking, 4 D x 4 times more efficient than a full MD simulation of 
a .D-dimensional cubic system. As we have already seen in one- and two-dimension cases, 
numerical results of our hybrid simulations show good agreements with those of CFD sim- 
ulations as far as Ai/i MD and Ax/Imd remain small, say At/t M D < 2 and Ax/Imd < 2. In 
fact, the normalized standard deviation, f L dx J T dt(v — vns) 2 /TL, of the velocity profiles 
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FIG. 7: Time evolutions of the velocity profile for the cavity flow with Re=980. The left col- 
umn shows the Hybrid simulations and the right column shows the corresponding CFD results. 
Simulation parameters are summarized as SP VII in Table HI 

of hybrid method, v, and those of CFDs, t>Ns ; are l ess than 0.02 and 0.07 for the cases of 
Fig. [3] (a) and for Fig. [3] (b). As the ratios increase, solutions of our hybrid model start 
to fluctuate around the corresponding CFD results. The deviation becomes about 0.6 in 
the case of Fig. H] (a). It is worth mentioning that the instantaneous velocity fluctuations 
are notable at each time step, however, they can be removed almost perfectly by taking 
time averages. In the following part, we will discuss the nature of the fluctuations more in 
detail to examine possibilities of effectively controlling them in our future simulations where 
correct thermal fluctuations will be included. 

To handle the statistical noise explicitly, we rewrite Eq. ffTTl) as 

p = -p/ + e^(T: + i?')0, (is) 
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TABLE I: Simulation parameters. 



where the off-diagonal stress tensor T', which is to be determined by MD sampling, is 
decomposed into the non-fluctuating stress T[ and the fluctuating random stress R' due to 
the thermal noise. The magnitude of each component of the random stress included in MD 
sampling (-Rmdpij)> where p and q represent the index in Cartesian coordinates and do not 
follow the summation convention, should depend both on the size of the MD cell Imd and the 
length of time £ M d over which average is taken at the MD level; (i?MDp 9 ) = (^W^MD, ^md) 2 ), 
where R(l, t) represents the random stress tensor averaged in a cubic with a side length / 
and over a time duration t. 

At the CFD level which is discretized with a mesh size Ax and a time-step At, the phys- 
ically correct magnitude should be (-RcFDpg) = (-Rpg(Aa;, At) 2 ). If the central limit theorem, 
(R pq (l,t) 2 )oc l/l D t is assumed, the following simple formula can be used. 




This finally leads to the following very useful expression for the correctly fluctuating stress 
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tensor P, 

P = -pI + 

to be used in CFD instead of Eq. (jlip . This equation indicates that if we can re- weight ran- 
domly fluctuating part B! while the non-fluctuating part T[ being untouched, hydrodynamic 
simulations including correct thermal fluctuations can be done for complex fluids within the 
present framework. 

We note that the important key toward the development of fluctuating hybrid simulation 
is the separation of T[ and R' . We thus carried out spectral analysis for the fluctuations in 
the total stress tensor computed directly from MD simulations T' = T[ + R '. The discrete 
Fourier transformation of T' xy is defined as 

1 2A/-12A/-1 

K y {k} = 4^ E E ZuW ex PH fc " x), (16) 

n x =0 n y =0 

where x = (n x Ax,n y Ax) is the position of each lattice node (n x ,n y ), k = 
(2irm x /L, 2itm y /L) is the wave vector, integers, M is the lattice num- 

ber in each x- and y-axis, and T' xy {x} is defined as T' xy {x}=T' xy (x + Ax/2,y + Ax/2) for 
< x,y < L, f' xy {x}=T' xy {2L - x, y} for L < x < 2L, and f xy {x}=T xy {x,2L - y} for 
L < y < 2L. 

The power spectra (\n' xy {k}\ 2 ) calculated from our hybrid simulations of driven cavity 
flows are plotted in Fig. [8] (a) for the case of At/^MD = Ax/Imd = 1- This corresponds to 
the case of Fig. [H] (a). The angle bracket (■ • • ) means the time average taken in the steady 
state at CFD level. One can see that the overall structure is rather simple. There exists 
a relatively large peak around k = and rather flat distributions throughout the k plane. 
The former corresponds to the contributions from the non-fluctuating part and the later 
corresponds to the contributions from the random stress R'. The same quantity obtained 
by conventional fluctuating hydrodynamics using a constant Newtonian viscosity and the 
random stress whose intensity is determined by the fluctuation-dissipation theorem^ is 
shown in Fig. [S] (b) for a comparison.— Those two plots are surprisingly similar to each other 
including the fluctuation part. This means that our hybrid simulation generates fluctuations 
quite consistent with the fluctuating hydrodynamics with fluctuation-dissipation theorem in 
the case of Ax//md=^/^md=1- 
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FIG. 8: The fluctuations of T' xy for the case of cavity flow with Re=59. The power spectra 
(|7Z^{fc}| 2 ) for the present multi-scale model with Ax/Imd = Ai/*MD = 1 is shown in (a) and the 
corresponding result from the fluctuating hydrodynamics is shown in (b) for a comparison. II' 
represents the discrete Fourier transform of T'. m a is defined as m a =(L/27r)k a , where k is the 



wave vector. The insets on each figure shows the {\II' \ 2 }-m x plane. 



-xy 
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Next, one see how the fluctuations depend on the ratios Aa;/I MD and Ai/t M D in Fig- 
[9j Here, comparing to the reference case (a) [Aa;/Z MD =At/tMD=2], the number of particles 
used in MD simulations are doubled in the case of (b) [Ax/Imd—^-59, A£/£md=2], and both 
the number of particles and the sampling duration to take time average are doubled in the 
case of (c) [Ao;/7md=1-59, At/tuo—i]- It is seen that the noise intensity decreases with 
decreasing ratios Ax/Imd and At/tMD in a consistent way to the central limiting theorem 
Eq. (|14p i.e., the noise intensity in (b) is about a half of that in (a), and the intensity in (c) 
is about one fourth of that in (a). 

Finally, we mention other recently proposed methods based on a similar idea. In the refer- 



ence 



11 



a hybrid method is proposed for bulk and boundary problems. Several problems for 
one- or two-dimensional flows of simple Lennard- Jones and dumb-bell liquids are considered. 
We note that the present multi-scale hybrid method is different from the methods proposed 
in those references particularly on the constructions of the stress tensor. In our method, a 
rotation matrix which effectively transforms the tensors in the Cartesian coordinates used in 
CFD and MD simulations. We also replace the isotropic part of the stress tensor calculated 
by MD simulations with the pressure imposed by the incompressible condition in CFD. More 
specifically, only the pure shear stress is passed from MD to CFD for numerical efficiency 
and consistency. 
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FIG. 9: The fluctuations of T' xy for the case of cavity flow with Re=235. The power spectra 
(\n' xy {k}\ 2 ) is plotted in (a) for the case of Fig. \6\ (b). Only the number of particles are doubled 
in (b), while other parameters are unchanged from (a). In (c), both the number of particles and 
the sampling time of T' are doubled. II' represents the discrete Fourier transform of T' xy . m a 
is defined as m a =(L/27r)k a , where k is the wave vector. The insets on each figure shows the 
(\n' xy \ 2 )-m x plane. 

V. SUMMARY 

We proposed a multi-scale method for hybrid simulations of MD and CFD. Our method 
is based on direct computations of the local stress by performing non-equilibrium MD sim- 
ulations according to the local flow field at all lattice nodes of CFD. The validity of the 
method is tested by comparing the numerical results obtained by our method and usual 
CFD. We found that the results obtained by our hybrid method agree well with those of 
usual CFDs with the Newtonian constitutive relation when the mesh size and the time-step 
of CFD are not too large comparing to the cell size and sampling time of MD simulations. 
When the ratios At/t MD and Ax/Zmd become large, there appear large fluctuations in flow 
field of our hybrid simulations. It was, however, clarified by the spectral analysis that the 
stress tensor T" computed by MD simulations has a very simple structure. It is composed of 
the non-fluctuating component T[ and the random component R' which seem to obey simple 
central limiting theorem according to the system size and the duration of the MD sampling. 
We confirmed that the power spectrum of the non-fluctuating component Tl in MD sampling 
agrees well with that computed in usual CFD without fluctuation. The power spectrum of 
R' also showed a good agreement with numerical results of the fluctuating hydrodynamics 
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which obeys the fluctuation-dissipation theorem. 
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